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AHHOTaUlHSI 

The percolation threshold and wrapping probability Rod for the two- 
dimensional problem of continuum percolation on the surface of a Klein 
bottle have been calculated by the Monte Carlo method with the Newman- 
Ziff algorithm for completely permeable disks. It has been shown that 
the percolation threshold of disks on the Klein bottle coincides with the 
percolation threshold of disks on the surface of a torus, indicating that this 
threshold is topologically invariant. The scaling exponents determining 
corrections to the wrapping probability and critical concentration owing 
to the finite-size effects are also topologically invariant. At the same time, 
the quantities Roo are different for percolation on the torus and Klein 
bottle and are apparently determined by the topology of the surface. 
Furthermore, the difference between the Roo values for the torus and Klein 
bottle means that at least one of the percolation clusters is degenerate. 
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1 Introduction. Formulation of the problem 

Problems of the percolation theory are widely used both in fundamental physics [1- 
21] and in applications. The fundamental interest in the percolation theory 
is associated with studies of the conformal field theory, probability theory, 
statistical physics, and theory of random graphs [1,2,4,7,8,10-13,18-21]. In 
applied problems, models of the percolation theory are used to describe granulated 
and composite materials [22-24], propagation of illnesses [25,26], and reliability 
of networks [18,19,27,28]. Models of the percolation theory are the main models 
for the description of processes in structurally disordered porous media filled 
with liquids and gases [29-31]. The states and properties of such systems have 
been actively studied in recent years [32-48] . Phenomena beyond the traditional 
notions were detected for such systems. In particular, the disordered system 
of pores of a porous medium filled with a nonwetting liquid can undergo a 
dispersion transition, when the nonwetting liquid transits to an effectively «wetting» 
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state at the variation of the temperature and degree of filling of the porous 
medium [45,46] . The state of such a system can be nonergodic with an anomalously 
slow relaxation of nonequilibrium states [47,48[. 

A physical reason for the dispersion transition is the appearance of the 
collective «multiparticle» interaction of liquid nanoclusters in neighboring pores 
with various sizes [47,48] owing to the formation of the ground state in the form 
of the percolation cluster of liquid-filled pores inside the percolation cluster of 
empty pores of the porous medium. The relaxation of the system is a process 
of its successive transition through the local maxima of the energy of local 
configurations of clusters of empty and filled pores and depends on the degeneracy 
of the ground state [47,48] , which is determined by the percolation cluster of 
filled pores. It is known that continuum percolation models [29-31,47-49] are 
used to describe the behavior of a liquid under the confinement conditions in 
disordered porous media. For this reason, the problem of the universality and 
properties of the percolation cluster in the case of continuum percolation arises, 
including the problem of the degeneracy of the percolation cluster in this case. 
We note that similar problems of the universality of various characteristics in 
two-dimensional percolation theory were previously studied both analytically 
and numerically [12-15] for a sphere [12], a torus [14], and a Mobius strip [13] . 

In this work, the problem of the possible degeneracy of the percolation cluster 
in the case of continuum percolation is solved by the numerical Monte Carlo 
simulation. To this end, we calculate the wrapping probability i?oo in the two- 
dimensional problem of percolation on a Klein bottle and compare with the 
results of the known calculation of percolation on the torus [14]. It is shown 
that the percolation thresholds and scaling exponents for the torus and Klein 
bottle coincide with each other, whereas the quantities it!co are significantly 
different. The last circumstance can be due to the degeneracy of the percolation 
cluster on one or both surfaces. 


2 Method of the calculation 

The Monte Carlo calculations were performed on square systems of various sizes 
with the Newman-Ziff algorithm [50] . This algorithm is currently the fastest 
algorithm of the search for the percolation threshold and ensures the calculation 
time proportional to the number of objects in a system (e.g., the calculation 
time in the Hoshen-Kopelman algorithm [51] is proportional to the square of 
the number of objects in the system). The main idea of the algorithm is that 
objects are added to the system one-by-one; after that, the conditions of the 
overlapping of an added object with already existing objects are verified, objects 
are joined into clusters, and, finally, percolation conditions are verified. The 
linear dependence on the number of objects is reached owing to the use of the 
union find with path compression algorithm [52,53] for the joining of clusters. 

In this work, completely permeable disks with the unit diameter d = 1 are 
used as objects. According to the Newman-Ziff algorithm, the following actions 
are performed in each calculation iteration: 

1. The addition of a disk to the system. 

2. The test of the conditions of overlapping of the added disk with the disks 
already existing in the system and joining of disks into clusters. 
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3. The test of the percolation conditions. 

The addition of a disk to the system is reduced to the generation of two 
random numbers determining the position of its center. The random numbers 
were generated according to the Mersenne-Twister algorithm having a period 
of — 1. 

The test of the conditions of overlapping of the added disk with the disks 
already existing in the system is complicated because the number of neighbors 
of the disk in the case of continuum percolation is unknown. In order to ensure 
the effective operation of the algorithm, we used an approach proposed in [14]; 
the essence of this approach is as follows. The calculation domain (a square 
with a certain size) is divided into squares with a side equal to the diameter of 
a disk, i.e., with a unit side. A disk located in a certain cell can overlap only 
with disks located in the same cell or in eight cells surrounding the given cell. 
Thus, the condition of overlapping is verified only for disks in the given cell and 
eight neighboring cells. 

Disks are joined into clusters according to the Newman-Ziff algorithm [50]. 
For the subsequent test of the percolation conditions, the shifts over two axes 
from the given disk to the root disk of the cluster are calculated in the process 
of joining. This path is formally a floating-point number, but it is easier to 
use the method proposed in [14] and to calculate the path from the square cell 
containing the disk to the cell of the root disk, rather than the path from the disk 
to the root of the cluster. This method allows using integers in the calculation 
of the path and in the test of the percolation conditions. 



Pnc. 1: Gluing of the boundaries providing the (left) torus and (right) Klein 
bottle. 

Particular attention should be paid to the formulation of boundary conditions, 
because they determine the topology of the surface on which the calculation is 
performed. The standard periodic boundary conditions correspond to percolation 
on the torus: the points of the boundaries in the vertical and horizontal directions 
are identified as is shown in Fig. 1 on the left. This corresponds to the «gluing» of 
a square sheet providing the torus. The gluing according to Fig. 1 on the right 
provides the Klein bottle. In terms of boundary conditions, this means that 
points on the upper and lower sides of the square are identified in the «mirror» 
manner, which creates additional difficulties for the test of the percolation 
conditions. 

In order to verify the percolation conditions in each of the directions, we 
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used the method proposed in [54] adapted for continuum percolation [14]; the 
essence of this method is as follows. Each addition of a new disk to the system 
is accompanied by the procedure of joining of clusters. If the added disk (disk 
1) overlaps with another disk (disk 2) and both disks belong to one cluster, 
there are two paths to their common root disk: through disk 1 and through disk 
2. If the sum of these paths in each of the directions is no more than unity, 
percolation does not occur. If the sum of these paths in a certain direction is 
larger than or equal to the size of the system in this direction, the wrapping 
cluster in this direction appears. 
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Phc. 2: Example of the square calculation domain for L = 4. Digits and blue 
color mark disks in the lower square, whereas primed digits and red color indicate 
their reflections in the upper square with allowance for mirror reflection in the 
vertical axis. 

However, it is noteworthy that the described method can be used to verify 
the percolation conditions only in the case of periodic boundary conditions and 
percolation on the torus. This method is inapplicable to the Klein bottle because 
gluing in one of the directions is mirror and the path to the root disk at the 
transition through the corresponding boundary can change by a value different 
from 0 and ±1. This difflculty is overcome as follows. Instead of the calculation 
in the square calculation domain with the size L, we performed the calculation 
in a rectangular calculation domain with a vertical size of 2L and a horizontal 
size of L, i.e., in a rectangle consisting of two squares in the vertical axis. When 
adding a new disk, two random numbers x, y were generated in the interval 
O^L, the disk was added to the lower square of the calculation domain, clusters 
were joined, and percolation conditions were verified. After that, a «copy» of 
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this disk was added to the upper square with allowance for mirror reflection 
in the vertical axis, i.e., with the coordinates L — x and y L\ then, clusters 
were joined, and the percolation conditions in the axes were verified. Figure 2 
exemplifies the used square calculation domain for the size L = A. Digits and 
blue color mark disks in the lower square, whereas primed digits and red color 
indicate their reflections in the upper square with allowance for mirror reflection 
in the vertical axis. 

The use of the rectangular calculation domain increases the calculation time 
with the algorithm as compared to the calculation in the square calculation 
domain, but eliminates the problem of verifying the percolation conditions, 
because the rectangular calculation domain has periodic boundary conditions in 
both axes and makes it possible to applied a method based on the summation 
of paths to the root disk. 

If each calculation ends at the number of disks n that results in the first 
appearance of the wrapping cluster, the probability of the existence of the 
wrapping cluster in the canonical ensemble (n) is equal to the fraction of 
the calculations ending at the number of disks no more than n. The wrapping 
probability Rl in the grand canonical ensemble is a function of the dimensionless 
concentration (degree of Ailing) 77 = nS/Li^ {S is the area of one disk) and is 
obtained by the convolution of the probability with the Poisson distribution 
with the mean value A = rjLp' /S [14]: 





( 1 ) 


At a fixed size of the system L, each calculation was performed until percolation 
occurred in both axes and the number of disks resulting in the first appearance 
of the wrapping cluster in each of the axes was fixed. Using these data and 
relation (1), we calculated six different probabilities Rl oi the appearance of 
the wrapping cluster: 

1. In any of the axes i?|. 

2. In both axes simultaneously 

3. In the horizontal axis R^ independently of the presence of a cluster in the 
vertical axis. 

4. In the vertical axis R^ independently of the presence of a cluster in 
horizontal axis. 

5. In the vertical axis, but not in the horizontal axis R]^. 

6 . In the horizontal axis, but not in the vertical axis R]^. 

It is worth noting that the vertical and horizontal axes are not equivalent in 
contrast to percolation on the torus. Consequently, in the general case, the 
probabilities are not equal to each other: ^ R^ and R^^ ^ At the same 

time, only three of six probabilities are independent because of the obvious 
relations 



( 2 ) 

(3) 
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Rf = R’l-Rl (4) 

The main problem in the simulation of percolation processes is the calculation 
of the percolation threshold rjc- The simulation of percolation on the Klein 
bottle is complicated because the asymptotic values of the probabilities i?oo 
are unknown in contrast to percolation on the torus. As a result, the relation 
Rl {vl) = Roc cannot be used to determine the dependence {L), which 
approaches rjc in the limit L —>■ oo. 

A method for the determination of the percolation threshold rjc based on 
a nonmonotonic behavior of the functions R\^ and R]^ was proposed in [50]. 
The positions of the maxima of the functions R}^ and R]^ should approach the 
percolation threshold rjc with an increase in the size of the system L ^ oo. Some 
other authors [55,56[ use the expression 

Rl [vl] = Rl/2 {vl) (5) 

to calculate the critical concentration at a given L value. The sequence of critical 
concentrations 77 ^ thus obtained should converge to the percolation threshold in 
the limit L —)■ c». In this work, we used both methods to estimate the critical 
concentration. 

To test the proposed simulation method, we calculated continuum percolation 
on the torus with the use of the reflection of the disk placed in the lower 
square to the upper square. In this calculation, in order to obtain the torus 
after gluing, reflection in the vertical axis of the disk in the upper square was 
not performed. The resulting critical concentration rjc and probabilities of the 
percolation transition R^o are very close to the results of work [14[. 


3 Results and discussion 

The calculations were performed for 21 sizes of systems in the range from L = 8 
to L — 3072. For each size, we performed N computer experiments: N > 10® 
for sizes 8 < L < 150, N > 10® for sizes 150 < L < 600, N > 10^ for sizes 
600 < L < 1200, and N > 10® for sizes L < 1200. Wrapping probabilities 
Rl iv) for some sizes from L = 8 to L = 3072 are plotted in Fig. 3. It is 
seen that R^ 7 ^ R'^ (Fig. 3a,b) and ^ R^ (Fig- 3e,f) because of the 
aforementioned nonequivalence of the vertical and horizontal axes in the case of 
the Klein bottle. 

Since the percolation threshold rjc for percolation on the Klein bottle is 
unknown, this threshold was calculated with the use of the dependences R^. In 
order to determine rjc, we obtained the dependence of the critical concentration 
on the size of the system rjL (A), which approaches the percolation threshold in 
the limit L —>■ 00 . The quantity (A) was determined by two methods: from 
the position of the maximum of R^^ and R^ and from relation (5). 

The rate of convergence ? 7 l to the percolation threshold rjc is determined 
by two factors: the characteristic width of the percolation transition region and 
the rate of convergence of the transition probability R^ to its asymptotic value 
Roc- The characteristic size of the transition region depends on the size of the 
system as , where v = 4/3 is the universal critical exponent for two- 

dimensional percolation. Our calculations show that, in the case of percolation 
on the Klein bottle, the characteristic size of the percolation transition region. 
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Phc. 3: Wrapping probabilities Rl for various sizes of systems from L = 8 to 
L = 3072: (a) R'l, (b) Rl, (c) Rl, (d) Rl, (e) , and (f) Rf. The dashed 

lines denote the expected values of rj^ and Rea- 


which is determined by the derivative Rj^ (vl), decreases really with an increase 
in L as This is seen in Fig. 4, where R^ ( 77 L) is plotted together with the 

straight line 0.361L^/'*. Three other derivatives R^i (rj^), {vl), and R"^ (rj^) 

have the same dependence. 

The rate of convergence of the probability of the percolation transition Rl 
to i?oo is determined by the boundary conditions [50] and for percolation on the 
torus depends on L as L~^. It will be shown below that the dependence L~^ [14] 
is also valid for percolation on the Klein bottle. Thus, the rate of convergence 
of the critical concentration r]L to the percolation threshold rjc in the case of 
percolation on the Klein bottle is determined by the dependence 

The semi-log plots of the critical concentrations tjl (L) derived from the 
maxima of R^} (open diamonds), R^ (closed diamonds) and from (5) for 
four probabilities: Rl (upward pointing triangles), R^ (circles), R^ (squares), 
and Rl (downward pointing triangles) are shown in Fig. 5. It is seen that 
all dependences with an increase in the size approach the same value. The 
percolation threshold obtained in this work is rjc = 1.128087(1), which almost 
coincides with the percolation threshold in the case of continuum percolation of 
permeable disks on the torus [14] . This indicates that the percolation threshold 
r]c is independent of the topology of the surface and is a topological invariant. 

The situation with the asymptotic probabilities of the percolation transition 
i?oo is different. These probabilities for percolation on the torus were calculated 
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Phc. 4: Derivative R\ (rj^) versus the size of the system L. 

in [57] from the conformal field theory and were confirmed in the numerical 
calculations [14,50[ for continuum and lattice percolation: 

R^= Rl^= 0.521058289..., 

= 0.690473724..., 

Rl^ = 0.351642853..., 

R^ = R^^ = 0.169415435.... 

These relations are invalid for percolation on the Klein bottle. To calculate the 
i?oo values, we used the following procedure. Taking into account the known 
percolation threshold rjc, we calculated the values Rl {rjc) and estimated the 
asymptotic value i?oo from this dependence by the least squares method with 
the function ^ + i?oo with the unknown parameters a, c, and i?oo- Among 
these parameters, the scaling coefficient c presenting the rate of convergence is 
of interest in addition to i?oo- It is noteworthy that the functions Rl in the 
case of percolation on the torus satisfy a power law with an exponent of-2; i.e., 
c = 2. The results of the calculations are summarized in the table 1. 

It is seen that the exponents c for all types of the probabilities are close to 2. 
This means that the scaling dependence of the probabilities i?L is independent of 
the topology, which confirms the above conclusion that the critical concentration 
rjL tends to the percolation threshold r]c according to the law This law 

is independent of the topology of the surface. The rate of convergence of the 
dependences Rl, which is determined by the coefficient c, is illustrated in Fig. 6, 



1.1292 



Phc. 5: Estimated critical filling factors for continuum percolation on a Klein 
bottle, derived from the maxima of (open diamonds), (closed diamonds) 
and from (5) for (upward pointing triangles), (circles), (squares), and 
(downward pointing triangles). Lines are 


Tadjiupa 1: Results of the calculations of the asymptotic values R^o by the least 
squares method with the function + R^o ■ The rightmost column presents the 
i?oo values for continuum percolation on the torus [50,57]. 


type 


c 

TDtorUS 

Roo 

TDh. 

0.48064(7) 

1.9(0) 

0.521058... 

T)V 

Roo 

0.51934(9) 

1.9(8) 

R%o 

0.68813(8) 

1 .8(8) 

0.690473... 


0.31185(8) 

1.8(3) 

0.351642... 

p7n 

0.16878(9) 

1.7(2) 

0.169415... 

Roo 

0.20749(1) 

1.8(9) 


which shows the dependences of the absolute value of the difference Rl — Roo 
on the size of the system. The errors on the plots were calculated from the 
relation [14] 


^Rl 


Rl iv) (1 - Rl jv)) 
N 


( 6 ) 


It is worth noting that the asymptotic values of the probabilities for percolation 
on the Klein bottle are smaller than those for percolation on the torus 1. 
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Phc. 6: Difference Rl — Roo versus the size of the system L for R^ (a), R" (b), 
i?® (c), R^ (d), (e), and R^^ (f). The solid lines are the dependences L~'^ 

with the exponents c from the table 1. 


Moreover, the probabilities and do not coincide with the probabilities 
R^ and because the vertical and horizontal axes are not equivalent in the 
case of the Klein bottle. We also note that the wrapping probability for the 
vertical axis on the Klein bottle R^ is closer to the wrapping probability on 
the torus than the wrapping probability for the horizontal axis R^. It should 
be reminded that the standard «gluing» without reflection corresponding to 
periodic boundary conditions is used for the horizontal axis, whereas «gluing» 
with reflection is applied for the vertical axis. At the same time, the probability 
is closer to the value on the torus than the probability R^. 

4 Conclusions 

The percolation threshold rjc and asymptotic values of the wrapping probabilities 
i?oo for the continuum percolation of permeable disks on the Klein bottle, as 
well as scaling exponents determining the contribution of finite-size effects to the 
critical concentration and wrapping probabilities Rl, have been obtained. 
It has been shown that the percolation threshold rjc and scaling exponents are 
independent of the topology of the surface and coincide with the respective 
values obtained for continuum percolation on the torus. At the same time, the 
wrapping probabilities differ from the respective values on the torus. In addition, 
since the vertical and horizontal axes are not equivalent in the case of percolation 
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on the Klein bottle, ^ and R^ ^ R^. In this case, the probability 
is closer to the value on the torus than the probability R!^, whereas the 
probability is farther from the probability on the torus than the probability 
R^. The difference between the probabilities Roo for the torus and Klein bottle 
means that the percolation cluster on one of these surfaces (or on both surfaces) 
is degenerate. 
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